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Nomenclature 


a speed of sound 

C v = 2(p ~^g o) surface pressure coefficient 

1 Poo U oo 

D cylinder diameter 

i =v / — T, imaginary number 

M Mach number 

p pressure 

Re = l V J D , Reynolds number 

u,v,w Cartesian fluid velocity components 
u* friction velocity 

V Cartesian velocity vector 

x,y,z Cartesian coordinates 
y+ = u*y_ nondimensional wall distance 

tf v ’ 


Greek: 

p fluid density 

0 azimuthal angle (x — y plane) 

v kinematic viscosity 

Superscript: 

' perturbation quantity (e.g. p' = p — p^) 

Subscript: 

oo freestream quantity 

rms root mean square 


Recent experience in the application of an optimized, second-order, backward-difference (BDF20PT) tem- 
poral scheme is reported. The primary focus of the work is on obtaining accurate solutions of the unsteady 
Reynolds-averaged Navier-Stokes equations over long periods of time for aerodynamic problems of interest. 
The baseline flow solver under consideration uses a particular BDF20PT temporal scheme with a dual-time- 
stepping algorithm for advancing the flow solutions in time. Numerical difficulties are encountered with this 
scheme when the flow code is run for a large number of time steps, a behavior not seen with the standard second- 
order, backward-difference, temporal scheme. Based on a stability analysis, slight modifications to the BDF20PT 
scheme are suggested. The performance and accuracy of this modified scheme is assessed by comparing the com- 
putational results with other numerical schemes and experimental data. 


I. Introduction 

Significant progress has been made in the field of computational fluid dynamics (CFD) during the last three decades, 
and modem CFD codes are applied routinely for solving steady, viscous flows on complex configurations using the Euler 
and Reynolds averaged Navier-Stokes equations. 1-6 However, the computational methodology is not as mature for the 
solution of unsteady flows, which is one of the next big challenges for the CFD community. The unsteadiness could be 
self-induced, such as for large separated flows, be created by flow control devices, or generated by many other mecha- 
nisms. Obtaining statistically meaningful solutions for such flow problems requires significant computational resources. 
However, advances in multi-core, multi-processor computer technology and numerical algorithms that take advantage of 
massive parallelization have resulted in significant reduction in wall clock time to obtain statistically converged solutions 
for problems of practical interest. 

An interesting and pertinent application of unsteady flow simulations is to provide input data to acoustic noise 
propagation codes. Such applications frequently have significant fluctuations over a range of frequencies and require 
unsteady solutions over a very large number of time steps to produce broadband statistical data with adequate temporal 
resolution. However, care must be exercised in selecting the numerical scheme for such applications because long time 
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computations allow minute numerical instabilities (e.g. unstable eigenvalues of magnitude 10 -5 ) to grow without bound, 
a phenomena that might remain unnoticed for short time computations. 

One possibility of improving the numerical efficiency is by adopting a higher-order, multi-stage, temporal scheme 
(e.g. 3rd- or 4th-order, implicit, Runge-Kutta 7 ). Indeed, the efficacy of high-order schemes has been demonstrated by 
several researchers on problems with varied error tolerances, particularly those that require high levels of accuracy. 8-10 
Adopting high-order schemes in large-scale engineering applications is not without complications. High-order schemes 
place an additional burden on the numerical algorithm used to solve the nonlinear set of equations, and only realize 
their full potential when the solver delivers a convergence rate that is relatively independent of time- step. Time- step 
independent convergence rates may be a difficult goal to achieve in general. Solver technology is one of the least mature 
disciplines among the many that are needed in the CFD community, and is presently the focus of extensive research. 
Placing any additional burden on the solver is not to be taken lightly. Hence, despite the recent successes of high-order 
schemes, there is still the need for reliable and robust, second-order methods. 

Herein, we focus strictly on multistep methods; in particular, the backward-differentiation formula (BDF) class of 
multistep schemes. These schemes have historically been very popular for large scale scientific computations because of 
their simplicity, nonlinear robustness, and relative efficiency. Perhaps the most popular of the BDF class of schemes is 
the second-order scheme that requires three time levels (BDF2). In addition to being second order accurate, it can be used 
for arbitrarily large time-steps (A-stability) without encountering numerical instability. Our ultimate goal in this work 
is to identify the most accurate BDF2 scheme that is stable for arbitrarily large time-steps. We identify an alternative 
second-order BDF2 scheme, (one of a class of optimized, second-order, backward difference methods referred to as 
“BDF20PT schemes”), with an error constant half as large as the conventional BDF2 scheme. The new scheme is 
extensively tested on model problems as well as a realistic problem involving the unsteady flow past a circular cylinder. 

II. Governing Equations and Numerical Scheme 

For the present work, the Reynolds-averaged Navier-Stokes equations are used as the governing equations to model 
the flow field. The discrete form of these equations on unstructured grids is solved using the Computational Fluid 
Dynamics (CFD) code known as FUN3D developed originally by Anderson and Bonhaus. 11 The FUN3D code has 
gone through significant modifications over the years by a team of researchers at NASA Langley using modem software 
practices. 12 Solutions are obtained either in a time-accurate manner with a constant time step at every grid point or with 
variable time stepping to accelerate convergence to a steady state. The effects of turbulence are modeled by using either 
the Spalart-Allmaras (SA) 13 model or with Menter’s two-equation shear stress transport (SST) 14 model. 

FUN3D has a wide range of flux schemes and flux limiters from which to choose. For the current work, Roe’s flux- 
difference splitting scheme 15 is used without a flux limiter. For steady flows, the solution at each time step is updated 
with a backwards Euler time-differencing scheme and the use of local time stepping. At each time step, the linear 
system of equations is relaxed in a red-black fashion with a point implicit procedure. 16 For unsteady flow problems, 
a dual-time- stepping scheme 17 is employed. A variety of time marching schemes are available in FUN3D, including 
a second-order, backward-differencing formulation (BDF2), and an optimized, second-order, backward-differencing 
formulation (BDF20PT). These will be discussed further in the next section. 

III. Temporal Discretization 

We focus strictly on multistep methods in this study; in particular, we examine the backward-differentiation for- 
mula (BDF) class of multistep schemes. These schemes, popularized by Gear, 18 are widely used for the integration of 
differential equations. The BDF schemes applied to the differential equation u t = f(u , £), can be expressed in the form 

k 

5j\7 3 U n +i — Atf n -\- 1* 

3=0 

In these equations, V 3 represents the backward difference operator and At represents the time step. The values of the 
coefficients Sj are easily derivable from a Taylor series expansion and are reported in Ref. 19 (pp. 365 — 366). 

It is well known that only the first- and second-order BDF schemes (BDF1 and BDF2) are A- stable, while BDF 
schemes of order three and higher are only A(a)-stable. Dahlquist’s famous theorem 20 states that 

Theorem 1 Let p represent the polynomial order of the temporal discretization. An A-stahle multistep method must be 
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of order p < 2. Furthermore, if the order is 2, then the error constant satisfies 


Here, C is the coefficient of the leading truncation error term. 

^-stability requires that all points in the left half of the complex plane (LHP), including the imaginary axis, are 
damped over one time step. A (a) -stability allows for growth in a wedge region of the included angle, (90 — a), that 
lies in the LHP between the imaginary axis and the line joining the origin and the tangency point of the root locus of the 
stability curve. More detailed discussion on A (a } -stability is available in Refs. 20 and 21. 

This second barrier of Dahlquist limits the utility of the BDF class of schemes for p > 3 as a general purpose 
scheme. A noteworthy exception is the case p = 3. For example, Melson et al. 22 showed that the BDF3 scheme was 
viable for many CFD applications. However, numerical instabilities were encountered by Vatsa and Carpenter 23 when 
using the BDF3 scheme to simulate the unsteady flow created by synthetic jets, thus underscoring the pitfalls of using a 
conditionally stable scheme. 

The stability characteristics for various BDF schemes are given in Table 1 . 


Method 

Levels 

Order 

A(a) 

BDF1 

2 

1 

A{ 90°) 

BDF2 

3 

2 

A(90°) 

BDF3 

4 

3 

A(86.03°) 

BDF4 

5 

4 

21(73.35°) 


Table 1. Properties of BDF methods. 


Given the “impossibility” of constructing A-stable multistep BDF methods higher than second order (of any number 
of levels), an obvious related question is “What is the most accurate A-stable BDF scheme that can be constructed?” 
Or more precisely, “Given an arbitrary number of time levels, what is the A-stable BDF scheme with a leading order 
truncation error closest to the C < — ^ barrier?” Asymptotically, one could imagine approaching the barrier as the 
number of time levels is increased. Clearly, large-scale computations place practical limitations (storage and work) on 
the required number of time-levels. Furthermore, preliminary investigation suggests that most benefits are achieved by 
including one or two additional time levels. Thus, in this work we focus exclusively on the four- and five-time-level, two 
parameter family of BDF schemes constructed from a linear combination of the BDF2, BDF3 and BDF4 schemes. 

We now begin by re-optimizing the BDF2 scheme subject to the A-stability constraint and four- time-levels. § A 
one parameter family of four-time-level, second-order, multistep schemes [BDF20PT(/5)j, is achieved by using a linear 
combination of the BDF2 and BDF3 schemes: 

BDF20PT(/5) = (/5) BDF3 + (1 - /5) BDF2 (1) 

Next, we adjust the parameter (3 to minimize the leading order truncation error of the BDF20PT(/5) scheme, subject 
to the constraint that the resulting scheme is A-stable. Note that adjusting the parameter /5 merely scales the coefficient 
of the leading order truncation term, C = fj, of the original BDF2 scheme by the factor (1 — /5). Thus, based on 
Dahlquist’s theorem, we should test the schemes on the interval 0 < /5 < |. 

The stability properties of multistep schemes are determined by assessing where in the complex plane the roots of 
the temporal polynomial lie. The root locus curves of the BDF class of schemes is given by 

k 

M = yAt 1 — exp (—iO)Y (2) 

j = i J 

For example, given the BDF2 scheme (k = 2), the real part of the root locus curve is defined by Re(p) = § — 
2cos($) + ^cos(26), which upon inspection reveals that Re(p) > 0 because the root locus curve never crosses the 
imaginary axis. Thus, the scheme is A- stable because all LHP eigenvalues are damped. 

Combining Eqns. (1) and (2) yields an equation for the root locus of the BDF20PT(/5) scheme: 

H | sin 3 4 (|) (4/5cos($) + 2/5 — 3) 

+ i\ sin(0)[2/5 cos(2$) — 3(2 p + 1) co s(6) +4/5 + 6] 

§Nyukhtikov et al. 24 presented a slightly different exercise when deriving their optimized second-order scheme, known as BDF20PT(/3 
Their primary focus was the accuracy of four-level schemes, and did not specifically require A-stability. 


(3) 

= 0.58). 
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and A-stability is guaranteed if Re(p) > 0 for all 0. Differentiating Re(p ) with respect to 0 in equation 3 and solving 
for all roots reveals that extrema exist at 

{x — ► mr; n — 0, ±1, ±2, • • •}, — ► ±2 cos -1 

The values of Re(p) at these extrema are 


zb 


2^ 


{ 


0 , 


8/3 (2/3 — l) 3 1 

3 ’ 12/3 2 j 


If all the extremum values are strictly non-negative, then the resulting scheme is ,4-stable. Thus, the BDF20PT(/3) 
family of schemes is ,4-stable for (3 on the interval 0 </?<§• 

Note that among the four- time-level A-stable BDF20PT(/3) schemes, the most accurate is achieved for the value 
(3 = |, which will be referred to as BDF20PT(4). For this case, the leading order truncation term is exactly half of 
the original BDF2 scheme: i.e. C = It is also noted that the BDF20PT scheme of Nyukhtikov et al. , 24 which was 
constructed using f3 = 0.58, is A(ct)-stable, with a = 89.89°, and hence is weakly unstable. 

The stability characteristics of the four- time-level BDF20PT( / $) family of schemes for two (3 parameters are listed 
in Table 2. 


Method 

Levels 

p 

A(a) 

BDF20PT(4) 

4 

0.50 

,4(90.00°) 

BDF20PT(/3 = 0.58) 

4 

0.58 

,4(89.89°) 


Table 2. Properties of four-time -lev el BDF20PT((3) methods. 

A similar exercise can be performed by parameterizing the BDF2 scheme with two additional time levels: 

BDF20PT(/?, 7 ) = 7 BDF4 + (3 BDF3 + (1 - (3 - 7) BDF2 (4) 

The A-stable scheme having the lowest error constant is achieved for the case 7 = 1 — and (3 = — | + 2\[2 which 
will be referred to by BDF20PT(5). 

Fig. 1 (a) shows five stability locus curves. The inner and outermost curves correspond to BDF2 and BDF4, respec- 
tively. It is observed from this figure that a small portion of the root locus of the BDF3 and BDF4 schemes lies in the 
LHP, which indicates the A(az)-stable nature of these schemes. The BDF20PT(4) and BDF20PT(5) schemes do not 
cross into the LHP, and are therefore indicative of an A-stable scheme, just as is the case for the BDF2 scheme. Fig. 1 
(b) shows an expanded view of the stability locus near the origin. Observe that except for the origin, the BDF20PT(4) 
scheme is bounded away from the imaginary axis while the BDF20PT(5) scheme touches the axis at precisely two 
additional points. 

This stability analysis is further substantiated with a numerical test. The model equation U t + U x = 0; 0 < x < 
1, t > 0 is integrated in time using the BDF20PT(/3) schemes for f3 = 0(BDF2), 0.45, 0.48, 0.50, 0.58, 1.0(BDF3) 
as well as the five level BDF20PT(5) scheme. The initial condition is a V (x, 0) = sin(27nr). The spatial derivative 
U x is discretized using in all cases a skew- symmetric, second-order, central-difference operator. The skew- symmetry of 
the matrix is maintained near the boundaries by adopting a periodic spatial domain, thereby guaranteeing that the semi- 
discrete eigenvalues are confined to the imaginary axis. This setting increases the probability that at least one eigenvalue 
will reside in the unstable lobes of the BDF20PT(/3) schemes for (3 > Fig. 2 shows a convergence study comparing 
the aforementioned BDF20PT(/3) schemes. The solution error is plotted as a function of the time-step (CFL = ^ with 
Sx fixed and St varying) on a logarithmic plot. All schemes converge over at least a portion of the time-step refinement, 
and the convergence rate is design order in all cases. The cases for which the parameter [3 exceeds the threshold (3 = \ 
are only conditionally stable, as is indicated by the sudden divergence of the solutions. It is also noted from Fig. 2, that 
the formal accuracy of BDF2 and BDF20PT(/3) schemes are the same (i.e. second order), whereas the actual error in 
the solutions for the BDF20PT(4) scheme is approximately half compared to the BDF2 solutions at a given time step. 

A question arises concerning the efficacy of the BDF scheme derived from the five-level case: BDF20PT(5). In this 
case, the stability locus curve touches the imaginary axis at two points in addition to the origin. Modes that coincide 
with these two points on the imaginary axis will be neutrally stable and could potentially lead to stability problems in the 
presence of algebraic perturbations resulting from inadequate convergence of the nonlinear system at each time-level. 
The stability characteristics of the BDF20PT(5) scheme remains the ongoing focus of investigation, as does the task of 
finding the general expression of the coefficients for an arbitrary number of time levels. 
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In summary, the BDF2 scheme has been re-optimized using slightly different criteria than those proposed by Nyukhtikov 
et al } 4 It is proven that the leading-order truncation error of the conventional BDF2 scheme can be decreased by a factor 
2 or - ~ 2.64 by combining it with the BDF3 and BDF3/BDF4 schemes, respectively. 

IV. Results 

For practical 3-D aerodynamic applications, it is important to select a temporal scheme that is efficient and robust, and 
does not adversely impact the performance of the base solver. Having examined various higher order temporal schemes, 
it was decided to use the BDF20PT scheme with four time levels described by Eqn. (1) in the previous section. This 
scheme requires the solution vector be stored at one additional time level compared with the standard BDF2 scheme; 
however, the lead truncation error term is nearly halved with only a nominal increase in the operation count. 

We now examine the turbulent flow past a circular cylinder at a free stream Mach number of 0.166 and at a Reynolds 
number of 166,000 based on cylinder diameter, D, to match the test conditions of the experimental study of Jenkins 
et al. 25 The turbulence model used for this case is based on the two-equation shear stress transport (SST) model of 
Menter. 14 It is well known that turbulence models developed for steady flows, such as the SA and SST models, are 
overly diffusive 26,27 and suppress the three-dimensionality of the flow in CFD simulations for unsteady flows that are 
accompanied by large reverse flow regions. The Detached Eddy Simulation (DES) model suggested by Spalart 28 has 
been used with reasonable success for the computation of unsteady flows past circular cylinders. 26,29-31 Others have 
used various forms of hybrid RANS/LES models for such computations. 32,33 In the present work, the SST model was 
modified such that the production terms were set to zero outside the boundary-layer region in the manner described by 
Khorrami et al. 21 and Lockard et al. 34 which results in a much less diffusive model and lower eddy viscosity levels in 
the wake regions. 

Wall temperature based on adiabatic flat plate conditions 11 is imposed on the surface of the cylinder along with the 
no slip conditions. Riemann-invariant based boundary conditions are imposed on the far-field boundary. In the span wise 
direction, periodic boundary conditions are imposed at 2 = 0 and 2 = to more closely simulate an infinite span 
problem. For the initial computations, the far-field boundary is set at a distance of 15 cylinder diameters away from 
the center of the circular domain and a span of 3 cylinder diameters is used. A structured grid consisting of 289 x 289 
points in the circumferential and radial direction was created in a two-dimensional plane, and extruded to produce 97 
equally spaced planes in the span wise direction to create the three-dimensional grid used in these computations. The 
grid spacing was clustered near the cylinder surface with a normal grid spacing equal to 0.0002 of a cylinder diameter 
at the surface which results in a value the turbulent parameter, y + of order 1 . A partial view of the grid in the end plane 
is shown in Fig. 3, where the grid clustering near the surface and wake region is clearly visible. The choice of grid 
clustering in the wake region was influenced by the Detached Eddy Simulation (DES) arguments proposed by Spalart 28 
and to replicate the grid topology suggested for DES simulations for cylinders by Travin et al. 26 and Vatsa and Singer. 31 

IV.A. Original BDF20PT(/3 = 0.58) scheme 

The first series of results reported here were obtained with the BDF20PT( / # = 0.58) scheme originally proposed by 
Nyukhtikov et al. 24 for improving the accuracy of BDF2 scheme. Although this scheme is A(a )~ stable and, therefore, 
conditionally stable, its stability range is much broader than the BDF3 scheme (see Tables 1-2), and it is more accurate 
than the conventional BDF2 scheme at nominally the same computational cost. This scheme was used originally by 
Nyukhtikov et al 24 for several turbomachinery problems and did not show any signs of instability. Vatsa and Carpenter 23 
have examined the performance of this scheme for the unsteady Navier- Stokes equations and demonstrated the improved 
accuracy of BDF20PT scheme compared to the BDF2 scheme for several test problems without encountering any 
numerical instability. The BDF20PT scheme has also been used successfully with the FUN3D code for several problems 
of aerodynamic interest by Biedron et al . 35,36 In the present work, we undertook the task of assessing the performance 
and accuracy of this scheme along with that of the standard BDF2 scheme for computing the unsteady flow past a circular 
cylinder over long periods of time with eventual application to acoustic noise prediction. 

The flow field was initialized as free stream and the FUN3D code was run for about 2000 iterations in a steady mode 
using the standard SST turbulence model. These solutions were specified as initial conditions for the unsteady flow 
computations performed with the modified SST turbulence model described earlier in the previous section. As can be 
seen from Fig. 4, care was taken to ensure that the residuals in the dual time-stepping converge 3-4 orders for the hy- 
drodynamic equations. The values in the plot are the L 2 norm of the subiteration residuals over the entire computational 
domain. The residuals for the turbulence equations only drop 2 orders of magnitude because the convergence is poor in 
highly unsteady regions away from the body where the modifications to the SST model minimizes the influence of the 
turbulence model. 
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Instantaneous contour plots of the density in the z=0 plane, based on flow solutions obtained after running 10K, 20K 
and 26K time steps (K represents a value of 1000 in this context, thus 10K stands for 10,000) with this BDF20PT(/3 = 
0.58) scheme are shown in Figs. 5(a)-(c), respectively. Fig. 5 (a) displays the typical vortex pattern behind the cylinder 
which is confined to the wake region and becomes weaker further away from the cylinder. The fluctuations in the flow 
field away from the cylinder and wake region are primarily acoustic and cannot be seen on this scale at this instant in 
time. However, when we examine the solutions at 20K time steps, the density contours indicate the formation of spurious 
waves. These waves get stronger with time (see Fig. 5 (c) at 26K steps), and eventually overwhelm the solution in the 
whole domain. 

Although the analysis in section III would suggest that the BDF20PT(/3 = 0.58) is admitting instabilities because it 
is only conditionally stable, the problem being solved here is much more complicated than the model problem used for 
the analysis. Hence, numerical experimentation was used to rule out other possible sources of the spurious oscillations. 
Several possible causes for the instabilities were investigated in a systematic fashion. Characteristic-based boundary 
conditions in the downstream region of the viscous wake can induce numerical instabilities, especially in subsonic 
flow regions when sufficiently fine grids are used near the downstream boundary. The combination of specifying the 
downstream pressure with several different forms of far-field characteristic conditions had minimal effect. The problem 
also persisted with tetrahedral grids, which are more common for unstructured grid codes. 

The distance of the far-field boundary was then doubled, i.e. it was placed at a distance of 30 cylinder diameters 
away from the center of the cylinder. We also coarsened the grid in the far-field to be more isotropic so as to reduce the 
reflections from the far field. A global view of this grid in the z=0 plane is presented in Fig. 6, where an inviscid type 
of grid is seen near the far-field boundary. The solutions obtained on this grid after 15K and 30K time steps are shown 
in Figs. 7 (a) and (b), respectively. Whereas the solution after 15K time steps is well behaved and displays the expected 
trends, the instabilities have polluted the flow structure completely by the later point in time at 30K time steps. 

Several other numerical studies were conducted to isolate the cause. We discovered that these instabilities could be 
suppressed by coarsening the spanwise grid distribution by a factor of 2. Conversely, if the spanwise extent was reduced 
by half, while still maintaining the number of points in the spanwise direction of the original grid, spurious oscillations 
were observed in the computational results after about 25 K time steps. Based on this exercise, we inferred that these 
numerical instabilities are more likely to appear for cases where the numerical dissipation is reduced via finer grids. 

As a final step to ensure that there were no other sources of instabilities in our numerical scheme, the computations 
were repeated with the standard BDF2 scheme, which is theoretically unconditionally stable or A-stable. Computational 
results after 50K time steps with the BDF2 scheme are shown in Figs. 8 (a)-(b), and these results display no sign of 
instabilities. 

Based on the results obtained in these tests, the BDF20PT( / # = 0.58) scheme appears to permit certain unstable 
modes to grow because of its conditional stability. One may not encounter such instabilities in the initial stages of a 
computation, or for problems where sufficient physical or numerical dissipation is present to damp these out. However, 
one should exercise caution when using this scheme for problems requiring integration over long periods of time. 

IV.B. Modified BDF20PT(/? = 0.48) scheme 

For the next series of solutions, we used the modified BDF20PT(/3) scheme with a value of [3 = 0.48. Based on the 
stability analysis presented in Section III of this paper, it is recognized that the BDF20PT schemes are A-stable for any 
value of (3 < 0.50. A value of 0.48 for the [3 parameter was chosen to be conservative. The resulting scheme is an 
A-stable scheme, and it reduces the temporal error by approximately a factor of 2 when compared with the standard 
BDF2 scheme. For remainder of this paper, solutions are presented on the grid where far-field boundary was placed at a 
distance of 30 D away from center (see Fig. 6). 

The instantaneous density contours in the z=0 plane after 50K time steps with the modified BDF20PT( / # = 0.48) 
scheme are shown in Fig. 9 (a)-(b). These results are free of any sign of instability and are qualitatively similar to the 
results obtained with the standard BDF2 scheme shown in Figs. 8 (a)-(b) for the same grid. We have run the code to 90K 
iterations in time to monitor the performance of this scheme, and have not detected any sign of numerical instability or 
spurious behavior in these solutions. 

The time-averaged pressure distributions on the cylinder surface were also examined. We started averaging the flow 
field data after running the code for 20K time steps after the initial transients had left the computational domain. The 
time-averaged surface pressure coefficient distribution at each circumferential location was also averaged in the spanwise 
direction (assuming spanwise homogeneity) to obtain the time-averaged data as a function of azimuthal angle, 6 . The 
angle is measured from the upstream stagnation point and is positive in the clockwise direction. The results from the 
BDF2 and BDF20PT(/3 = 0.48) schemes are compared with experimentally measured data 25,37 for the cylinder in Fig. 
10. Note that the flow was tripped using a transition strip between 6^ = 50 — 60° in the experimental case to simulate 
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turbulent transition on the cylinder surface, whereas the CFD code does not have the mechanism to force transition in 
this manner. Therefore, the CFD code was run in a fully turbulent mode. However, the actual transition in the CFD 
solutions occurred later compared to the transition trip location in the experiment. However, separation did occur before 
6 = 90° so that the flow became turbulent before separation. 

The pressure distributions from the two temporal schemes in FUN3D are almost indistinguishable in Fig. 10. The 
computational results compare well with the measured data near the stagnation region as well as on the aft side of the 
cylinder. However, the computational C p indicates an earlier separation, thereby resulting in under-expansion of the 
flow near the top and bottom of the cylinder. We believe that at this moderate value of Reynolds number, the CFD code 
with the SST turbulence model does not simulate the same transition to turbulent flow as in the experiment. Hence, the 
discrepancy in the pressure comparisons. We also obtained solutions for this case on the same grid with the structured 
grid flow solver CFL3D 38,39 using the BDF2 scheme, and the resulting pressure distributions from these computations 
are also shown in Fig. 10. The pressure distributions obtained from the CFL3D and FUN3D codes are in excellent 
agreement with each other. 

The distribution of the root-mean- square (rms) perturbed surface pressure coefficient ( Cp rrns ) is shown in Fig. 11. 
The FUN3D solutions obtained with the BDF2 and BDF20PT schemes are compared with the experimental data ob- 
tained in the Quiet Flow Facility at NASA Langley Research Center. 37 The unsteady pressure transducers on the cylinder 
were located every 45°. To obtained refined azimuthal resolution, the cylinders were rotated, and the transition strip had 
to be re-applied after each rotation. Such an arrangement can easily create some non-uniformity in the configuration and 
flow characteristics and introduce scatter in the measured data, as seen in Fig. 11. Again, the FUN3D results with the 
two temporal schemes compare well with each other and are in good agreement with the CFL3D results. The numerical 
solutions miss the peak in the rms because the flow separates too early in numerical simulations. However, most of the 
computed rms pressure distribution lies within the scatter of the experimental data. 

Comparing the efficiency of temporal schemes in real-world settings is often a very difficult task. For example, 
see Ref. 9 for a thorough discussion on this topic. To actually demonstrate any accuracy differences between the time 
integration methods for the cylinder shedding problem, a systematic time- step refinement would be required. However, 
the case is dominated by fluctuations at the vortex shedding frequency. Any improvement in the time integration would 
only have a minor effect on the time-averaged quantities examined here unless the time step was large enough to alter 
the shedding process. Theoretically, one should be able to use a larger time step with the improved scheme, but exactly 
how much larger is difficult to estimate because of the interaction of the temporal and spatial operators as well as the 
complicated equation sets that are solved in real-world problems. Most CFD codes are limited by spatial rather than 
temporal accuracy, so the optimal time step would be difficult to determine without a systematic refinement in space 
and time. However, such studies are often prohibitively expensive for unsteady flow problems. The temporal and spatial 
resolutions chosen are often based on formulas derived from simplified equation sets and engineering judgment. Hence, 
an improvement in the temporal accuracy is often used like a factor of safety. In this case, one is trading a little higher 
operation count and the storage of one additional solution vector for the improved factor of safety. 


V. Concluding Remarks 

A new class of optimized, second-order, backward-difference (BDF20PT) schemes has been developed that possess 
formal stability properties with only a slight reduction in the theoretical accuracy of previous optimized schemes. Sev- 
eral numerical studies were conducted to demonstrate that the BDF20PT(/3 = 0.58) scheme of Nyukhtikov et al . 24 does 
admit spurious solutions for long-time integration of the flow over a circular cylinder. This scheme is only conditionally 
stable and, under certain circumstances, permits the growth of instability modes. A slight modification of the scheme re- 
sults in an unconditionally stable scheme. Results obtained with the modified scheme were shown to be well behaved and 
free of any spurious oscillations. The proposed BDF20PT(/3 = 0.48) scheme still offers improved temporal accuracy 
compared with the standard BDF2 scheme at nominally the same computational cost, but with a higher storage require- 
ment (a 3 to 10% increase in the total memory requirement for typical CFD codes). Time-averaged solutions obtained for 
flow over an isolated cylinder were shown to be essentially identical between the BDF2 and new BDF20PT( / # = 0.48) 
schemes. Although the new temporal scheme is more accurate, to achieve a decreased computational cost through a 
larger time step, the interactions between the spatial and temporal operators must be evaluated in the context of the 
equation set and geometry being solved. Such studies are often prohibitively expensive to perform computationally, and 
analysis techniques are not general enough to handle most real-world problems. Not only are higher- accuracy temporal 
and spatial schemes needed, but also robust, accurate, and efficient error estimates. Only with improved error estimates 
can the full potential of higher- accuracy methods be realized. 
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(a) Full domain view 


(b) Expanded view near origin 


Figure 1. Stability locus for the BDF2 schemes 
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Figure 2. Convergence behavior of BDF2, BDF3 and BDF20PT schemes for test problem 
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Figure 3. Partial view of grid in z=0 plane 
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Figure 4. Convergence histories for dual-time stepping scheme 
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(b) 20K time steps 
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(c) 26K time steps 

Figure 5. Instantaneous density contours in z=0 plane at different times. Solutions obtained using the BDF20PT(/3 = 0.58) scheme 
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Figure 6. Global view of grid with far-field boundary at 15D in z=0 plane 




(a) 15K time steps (b) 30K time steps 

Figure 7. Instantaneous density contours at different times, far-field boundary at 30D. Solutions obtained using the BDF20PT(/3 = 0.58) 
scheme 
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(a) Global view (b) Expanded view near origin 

Figure 8. Instantaneous density contours at 50K time steps, far-field boundary at 30D. Solutions obtained using the BDF2 scheme 





(a) Global view (b) Expanded view near origin 

Figure 9. Instantaneous density contours at 50K time steps, far-field boundary at 30D. Solutions obtained using the BDF20PT(/3 = 0.48) 
scheme 


14 of 15 


American Institute of Aeronautics and Astronautics Paper 2010-0122 







Figure 10. Comparison of time-averaged surface pressure 



Figure 11. Comparison of perturbation pressure coefficient 
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